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The Michaelis-Menten equation has played a central role in our understanding of biochemical 
processes. It has long been understood how this equation approximates the dynamics of 
irreversible enzymatic reactions. However, a similar approximation in the case of networks, 
where the product of one reaction can act as an enzyme in another, has not been fully 
developed. Here we rigorously derive such an approximation in a class of coupled enzymatic 
networks where the individual interactions are of Michaelis-Menten type. We show that 
the sufficient conditions for the validity of the total quasi steady state assumption (tQSSA), 
obtained in a single protein case by Borghans, de Boer and Segel can be extended to sufficient 
conditions for the validity of the tQSSA in a large class of enzymatic networks. Secondly, we 
derive reduced equations that approximate the network's dynamics and involve only protein 
concentrations. This significantly reduces the number of equations necessary to model such 
systems. We prove the validity of this approximation using geometric singular perturbation 
theory and results about matrix differentiation. The ideas used in deriving the approximating 
equations are quite general, and can be used to systematize other model reductions. 

Keywords: Michaelis Menten, quasi steady state, total quasi steady state, protein 
interaction networks, coupled enzymatic networks, geometric singular perturbation 



1. Introduction 

The Michaelis-Menten (MM) scheme [3j [20] is a fundamental building block of many 
models of protein interactions: An enzyme, E, reacts with a protein, X, resulting in an 
intermediate complex, C. In turn, this complex can break down into a product, X p , and the 
enzyme E. It is frequently assumed that formation of C is reversible while its breakup is 
not. The process is represented by the following sequence of reactions [31 [201 EI] 



Frequently catalytic activity in protein interaction network is modeled by MM equa- 
tions [HJ EH E31 E2J El IH E51 E]- This gives rise to coupled enzymatic networks, where the 
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substrate of one reaction acts as enzyme in another reaction. A direct application of the 
law of mass action to such models typically leads to high dimensional differential equations 
which are often stiff, and difficult to study directly. 

A number of methods have been introduced to address these problems. Most of these 
methods are based on quasi steady state assumptions which take advantage of the differences 
in characteristic timescales of the quantities being modeled. It is typically assumed that the 
chemical species, or some combinations of chemical species, can be divided into two classes: 
One which equilibrates rapidly, and a second which evolves more slowly [121 US] • Assuming 
that the members of the first class equilibrate instantaneously leads to a reduced model 
involving only elements of the second class. 

Reduction methods differ in their assumption on which chemical species, or combi- 
nations thereof, are assigned to the two different classes. For instance, the standard quasi 
steady state assumption (sQSSA) posits that the concentrations of the intermediate com- 
plexes change quickly compared to the protein concentration pm EH EU El [23] . An alterna- 
tive is the reverse quasi steady state assumption (rQSSA) where the protein concentration 
is assumed to change rapidly compared to intermediate complexes [29]. Rigorous justifica- 
tions of these methods are largely available only for isolated reactions of the type shown in 
scheme (1.1), and the Goldbeter-Koshland switch^] [ID] . 



The total quasi steady state assumption (tQSSA) was introduced to broaden the range 
of parameters over which a quasi steady state assumption is valid. Under this assumption 
the concentration of the intermediate complex, C, evolves quickly compared to the sum 
of the intermediate complex and the protein concentration [21 [3H [181 ES [35] • Numerical 
experiments and heuristic arguments suggest that tQSSA may be valid in coupled enzymatic 
networks over a very broad set of parameters [5] . 

Here we aim to provide a theoretical foundation for the reductions used in numerical 
studies of enzymatic networks. A standard model reduction technique for systems involv- 
ing quantities that change on different timescales is geometric singular perturbation theory 
(GSPT) [7J, [HI H2]. For instance, this theory has been used by Khoo and Hegland to prove 
several results obtained earlier by Borghans, et al. using self consistency arguments [151 12]. 
GSPT has also been used to reduce other models of biochemical reactions [371 HH EE]- We 
derive a sufficient condition for the validity of tQSSA in arbitrary networks of proteins and 
enzymes provided the interactions are of MM type and can be modeled by mass action ki- 
netics. This directly extends previous work, like that of Pedersen, et al. [27] who proposed 
a sufficient condition for the validity of tQSSA in the Goldbeter-Koshland switch. 

The direct application of the tQSSA to coupled enzymatic networks generally leads 
to a differential-algebraic system. The algebraic part of this system consists of coupled 
quadratic equations that are typically impossible to solve. Our second aim is to show that, 
under certain assumptions on the structure of the network, it is possible to circumvent this 
problem using ideas introduced by Bennett, et al. [I]. This allows us to obtain reduced 



1 A Goldbeter-Koshland switch consist of two coupled reactions. One of these reactions frequently repre- 
sents protein phosphorylation, and the second dephosphorylation. 
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set of differential equations for a class of protein interaction networks in terms of protein 
concentrations only. 

We proceed as follows: In section [2] we review the original Michaelis-Menten scheme. 
We introduce terminology, and illustrate our approach in a simple setting. In this section 
we also give a brief overview of the theory of geometric singular perturbation theory, which 
is fundamental in proving the validity of the reduced equations. In section [3] we extend our 
approach to a well studied two protein network that plays part in the G2-to-mitosis phase 
(G2/M) transition in the eukaryotic cell cycle. We present the ideas in the most general 
setting in section [4], where we derive the general form of the reduced equations. Each section 
begins by the discussion of the tQSSA in the context of the network under consideration, 
and closes with a derivation of the reduced equations under the tQSSA, as well as sufficient 
conditions under which the tQSSA holds. A number of technical details used in the proofs 
of the main results are given in the appendices. We note that throughout the presentation 
the law of mass action is assumed to hold. 



2. Isolated Michaelis-Menten reaction 



The MM scheme is frequently used to model enzymatic processes in solution which are 
ubiquitous in biology. As discussed in the introduction, a number of different approaches 
have been proposed to justify the reduced equations mathematically. We start by giving a 
detailed overview of the tQSSA approach based on geometric singular perturbation theory 
(GSPT)[7]. The setting of a single MM type reaction will be used to introduce the main 
ideas and difficulties of reducing equations that describe larger reaction networks. 

For notational convenience we will use variable names to denote both a chemical species 
and its concentration. For instance, E denotes both an enzyme and its concentration. Re- 



action (1.1) reaction obeys two natural constrains: The total amount of protein and enzyme 
remain constant. Therefore, 



X + C + X p = X T , and E + C = E 7 



(2.1; 



for positive constants Xt and Et- In conjunction with the constraints (2.1), the following 



system of ordinary differential equations can be used to model reaction (1.1) 

-hX(E T -C) + fe_iC, X(0) = X T , 

faX(E T - C) - (jfe_i + k 2 )C, C(0) = 0. 



dX 

~dt 
dC 

~dt 



(2.2) 



2.1. The total quasi steady state assumption (tQSSA) 

Under the standard quasi steady state assumption (sQSSA), the concentration of the 



substrate-bound enzyme, C, equilibrates quickly, which allows system (2.2) to be reduced 



by one dimension. Sufficient conditions under which the sQSSA is valid have been studied 
extensively [10J, [30l E] . However, it has also been observed that the sQSSA is too restrictive [21 

M. 
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To obtain a reduction that is valid for a wider range of parameters, define X :— X + C. 



Eq. (\2.2\j can then be rewritten as 
-k 2 C, 

h[XE T — (X + Et + k m )C + C 2 ] 



dX 

~dt 
dC 

~dt 



X(0) 
(7(0) 



0. 



(2.3a) 
(2.3b) 



where k m = (k_i + k 2 )/ki is the Michaelis-Menten constant. 

The tQSSA posits that C equilibrates quickly compared to X [21 EI]. Under this 
assumption we obtain the following differential-algebraic system 



dX 
~dt 



-k 2 C, 



X{0) = X T , 



= h[XE T -(X + E T + k m )C + C 2 



(2.4a) 
(2.4b) 



Solving Eq. (2.4b) and noting that only the negative branch of solutions is stable, we can 



express C in terms of X to obtain a closed, first order differential equation for X, 
(X + E T + k m ) - y/{X + E T + k m f - AXE T 



dX 
~dt 



X(0) = x 7 



(2.5) 



Although the reduced equation is given in the X, C coordinates, it is easy to revert to the 



original variables X,C. Therefore, from Eq. (2.5) one can recover an approximation to the 



solution of Eq. (2.2) 



2.2. Extension of the tQSSA 



An essential step in the tQSSA reduction is the solution of the quadratic equation (2.4b) 



A direct extension of this approach to networks of chemical reaction typically leads to cou- 
pled system of quadratic equations [5j ETJ [26] . The solution of this system may not be 
unique, and generally needs to be obtained numerically. However, an approach introduced 
by Bennett, et al. [T], can be used to obtain the desired solution from a system of linear 
equations. 

In particular, we keep the tQSSA, but look for a reduced equation in the original 



coordinates, X, C. Using X = X + C to eliminate X from Eq. (2.4b), we obtain 

= h (X(E T -C)- k m C) . 



(2.6) 



Eq. (2.6) and Eq. (2.4b) are equivalent, but Eq. ( 2.6[ ) is linear in C, and leads to 



C 



XEj* 
k m ~\~ X 



and X = X + 



XEj- 
k m + X' 



Using these expressions formulas in Eq. (2.4a), and applying the chain rule gives 
XEt 



_d_ 
dX 



X + 



k m "I" X 



dX 
~dt 



XE- 



r 



k m + X 



dX 
~dt 



k 2 1 + 



k m Er 

(k m + xy 



XEt 

k m + X 
(2.7a) 
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The reduced Eq. (2.7a) was obtained under the assumption that there is no significant 
change in X = X+C during the rapid equilibration . Aft er equilibration, C = XEj-/(k m +X) 
(See Fig. [l]). Therefore, the initial value for Eq. (2.7a), denote by X(0), can be obtained 
from the initial values X(0),C(0) using 



X(0) + 



E T X (0) 
X (0) + k n 



X(0) + C(0) = x T . 



(2.7b) 



Fig. ^p) shows that the solutions of the full system (2.2) and the reduced system (2.7a) 
are close when initial conditions are mapped correctly. 
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Figure 1: Proper choice of the initial values of the reduced system. The empty circle at X = 1, C = 0, 
represents the initial value for the full system. The solid dot is the initial value of the reduced system. The 
dash-dotted (red) line represents the attracting slow manifold, (a) The solid curve represents the numerical 
solution of Eq. (2.3 1. The solution rapidly converges to the manifold, and evolves slowly along the manifold 
after this transient. The dashed line satisfies X = Xt- The solid dot at the intersection of the dashed 
line and the slow manifold represents the projection of the initial condition onto the slow manifold given by 
Eq. (2.4b). Thus X(0) = Xt is the proper initial condition for the reduced system ( |2.5[ ). (b) The solid line 
represents the numerical solution of Eq. (2.2). After a quick transient, the solution again converges to the 
slow manifold. However, since the initial transient is not orthogonal to the X axis, the initial conditions do 
not project vertically onto the slow manifold. Instead, the initial transient follows the line X + C — Xt 
(dashed) , and the intersection of this line and the slow manifold represents the proper choice of the initial 
value for Eq. (2.7a). (c) Comparison of solutions of Eq. (2.2) and the reduced system (2.7a I. The graph in 
the inset offers a magnified view of the boxed region, showing the quick transient to the slow manifold. We 
used: Xt 



E T = h 
system, X(0) = 0.83. 



&2 = 1) k-i = 3, which, using Eq. (2.7b), gives the initial condition for the reduced 
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The tQSSA implies that Eq. (2.3) can be approximated by Eq. (2.4). Therefore, to 
explore the conditions under which Eq. (2.7a) is a valid reduction of Eq. (2.2) we need 
to provide the asymptotic limits under which the transition from Eq. (2.3) to Eq. (2.4) 
is justified. Different sufficient conditions for the tQSSA have been obtained using self- 
consistency arguments [21 [M]. We follow the ideas of pairwise balance to look for a proper 
non-dimensionalisation of variables [3T1 [8] . Although this method gives a weaker result than 
the one obtained in |34|, it is easier to extend to networks of reactions. 



2.3. Review of Geometric singular perturbation (GSPT) 

Since geometric singular perturbation theory (GSPT) is essential in our reduction of 
the equations describing coupled enzymatic reactions, we here provide a very brief overview 
of the theory. Further details can be found in |36l [161 El EE] • Readers familiar with GSPT 



can skip to section 2.4 



Consider a system of ordinary differential equation of the form 

du 
e— 
dt 

dv 

v 



dt 

l k and v G I 



f(u,v,e), 
g(u,v,e), 



u(0) 
«(0) 

v e 



UQ, 



(2.8a) 



where u G M fe and u G K' with k, I > 1, and Uq G Vq G M! are initial values. The 
parameter e is assumed to be small and positive (0 < e < 1), the functions / and g smooth, 

0. (2.8b) 

The variable u is termed the fast variable, and v the slow variable. 

Assume that M := {(u, v) G R k+l | f(u, v, 0) = 0} is a compact, smooth manifold with 



f(u, v , 0) ^ 0, g(u, v, 0) ^ 0, and limeg(u,v,e) 



9/, 



U,V, 



inflowing boundary. Suppose further that the eigenvalues Aj of the Jacobian 
all satisfy Re(Xi) < 0, so that A4q is normally hyperbolic. Then, for e sufficiently small, the 

v an initi 

f(u,v,0) 



u(0) = u , 

0, v(0) = Vq, (2.9) 

es. After this transient, the solutions are (9(e) close to the solutions of the reduced 



solutions of Eq. (2.8a) follow an initial transient, which can be approximated by 

du 

ds 
dv 

ds 

where t 
system 

= f(u,v,0), 
dv 

— = g{u,v,0), v(0)=v o . (2.10) 

More precisely there is an invariant, slow manifold Ai e , 0(e) close to .Mo- Solutions of 
Eq. (2.8a) are attracted to Aio exponentially fast, and can be approximated by concatenating 
the fast transient described by Eq. (2.9), and the solution of the reduced Eq. (2.10). 



The slow manifold, Aio, consists of the fixed points of Eq. (2.9). The condition that 
the eigenvalues, Aj, of the Jacobian t£(u,v,0)\m o all satisfy Re(Xi) < implies that these 
fixed points are stable. 
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2.4. Validity of the tQSSA 



We next show that GSPT can be applied to Eq. (2.3), after a suitable rescaling of 
variables [3TJ [S] • Let 

t X(t) C(t) 



x{t) 



Xn 



c r 



(2.11) 



We have some freedom in defining /3 and T%- Using the method of pairwise balance [HI EI], 
we let 

XtEt Xt 



(3 



X T + E T + k r , 



and T 



x 



k 2 /3' 



In the rescaled variables, Eq. (2.3) takes the form 

dx 



Ei 



dr 
dc 



h (E T + X T + k m ) 2 dr 
Define the parameter 



-c, 



x 



X T x + E T + k Ti 
X T + E T + k m 



-c + 



XtE' 



rpJ2jrp 



ko 



(E T + X T + k % 



Eq 



x{0) 
c(0) 



k x T x X T E T k x (E T + X T + k Tl 



(2.12) 

1, (2.13a) 
0. (2.13b) 

(2.14) 



For small e, Eq. (2.13) is singularly perturbed and has the form given in Eq. (2.8a). Indeed 



we can apply GSPT to Eq. (2.13) directly since in the limit e — > the right hand side of 



Eq. (2.13) remains 0(1). Indeed, the requirement < e <C 1, is equivalent to the sufficient 
condition for the validity of the tQSSA derived in [2]. 



GSPT implies that for small e, solutions of Eq. (2.13) are close to those of the reduced 
system 

dx 
dr 



-c, 



x(0) 







x 



X T X + E T + kr, 
X T + E T + k m 



■c + 



XtE' 



T-LLiT 



(E T + X T + k r , 



(2.15a) 
(2.15b) 



The normal hyperbolicity and stability of the manifold defined by Eq. (2.15b) can be 



verified directly, and also follow from the results of section |4| It follows that GSPT can be 
applied to conclude that GSPT implies that Eq. (2.15) is a reduction of Eq. (2.13). 



The validity of the reduction in these rescaled equations implies its validity in the 



original coordinates: Eq. (2.13) is equivalent to Eq. (2.3) via the scaling given in Eq. (2.11). 



Hence, Eq. (2.4) and Eq. (2.15) are related by the same scaling relationship. We note that 



choosing the initial values of the intermediate complexes to be zero, implies that solutions 



of (2.13) remain 0(1) for small e (see section 4.6 for a detailed discussion). It follows that 



Eq. (2.4) is a valid reduction of Eq. ( |2.3 ) when e is sufficiently small. Hence, for e in the same 
range, Eq. ( |2.7a[ ), with initial values satisfying Eq. (2.7b), is a valid reduction of Eq. (2.2). 



Lemma [7] in the Appendix shows that e is always smaller than 1/4. Although this is 
suggestive, GSPT only guarantees the validity of the reduced equations in some unspecified 
range of e values. 
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3. Analysis of a two protein network 

We next show how the reduction described in the previous section extends to a network 
of MM reactions. Here the substrate of one reaction acts as an enzyme in another reaction. 
To illustrate the main ideas used in reducing the corresponding equations, we start with a 
concrete example of two interacting proteins. 
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(a) 



(b) 



® 



© 

A 



© 




Original system 
Reduced system 



time (t) 



Figure 2: A simplified description of interactions between two regulators of the G2-to-mitosis phase (G2/M) 
transition in the eukaryotic cell cycle [25] (See text), (a) X and Y phosphorylate and deactivate each 
other. For instance, the protein X exists in a phosphorylated X p and unphosphorylated X state, and the 
conversion X to X p is catalyzed by Y p . The con vers ion of X p to X is catalyzed by the phosphatase E\. 
(b) Comparison of the numerical solution of Eq. (3.1| and Eq. ( |3.8[). Here k% = 5, = 1, &2 = ^,E^ = 
10, El = 2,X T = 10, Y T = 10.1 as in [5 . The initial values for Eq. (|3.1|) are X(0) = 10, Y(0) = 1.1, X p (0) = 
0,r p (0) = 9,C X (0) = 0,C y (0) = 0,C|(0) = 0,C«(0) = 0,£i(0) = 10,^2(0) = 2. The initial values of the 
reduced system, X P (Q) = 0.12, Y p (0) = 0.83 are obtained by the projection onto the slow manifold defined 
byEq. 



Coupled enzymatic networks - tQSSA 



11 



Fig. |2p,) is a simplified depiction of the interactions between two regulators of the G2- 
to-mitosis phase (G2/M) transition in the eukaryotic cell cycle [25]. Here, Y represents MPF 
(M-phase promoting factor, a dimer of Cdc2 and cyclin B) and X represents Weel (a kinase 
that phosphorylates and deactivates Cdc2). The proteins exist in a phosphorylated state, 
X p , Y p , and an unphosphorylated state, X, Y, with the phosphorylated state being less active. 
The proteins X and Y deactivate each other, and hence act as antagonists. In this network 
Ei and E 2 represent phosphatases that catalyze the conversion of X p and Yj, to X and Y, 
respectively. Each dotted arrow in Fig. |2^i) is associated with exactly one MM type reaction 
in the list of reactions given below. The sources of the arrows act as enzymes. Therefore, 
Fig. [2^,) represents the following network of reactions 

kl k 2 kl e k 2 

Y p + X C x — > X p + Y p , Ei + X p ^ C x — y X + Ex, 

X p + Y hc y ^Yp + X p , E 2 + Y p hc e y ^Y + E 2 . 
k-i k-i 

To simplify the exposition, we have assumed some homogeneity in the rates. Since the total 
concentration of proteins and enzymes is assumed fixed, the system obeys the following set 
of constraints 

X T = X(t) + X p (t) + C x (t) + C y (t) + C e x {t), El = C e x (t) + Ei(t), 

Y T = Y(t) + Y p (t) + C x (t) + C y {t) + C c y (t), E T 2 = C e y (t) + E 2 (t), 

where Xt,Yt, El , El are constant and represent the total concentrations of the respective 
proteins and enzymes. Along with these constraints the concentrations of the ten species in 
the reaction evolve according to 



dX p 
~dT 



-ki (Yp -Y p - C x -C v - C e y ) Xp - hX p (El - C e x ) +k.iC e x + (*_! + k 2 )C y + k 2 C x , 



--Y =Ei 



dY 

= -ki(X T -X p -C x -C y - Cl)Yp - hY p (El - C e y )+k_iC e y + (k-i + k 2 )C x + k 2 Q 

dC x 



-v- 



dt 



dCy 

~dT 



=X =E 2 

h (X T -X p - C x -C y - CI) Y v - (*_! + k 2 )C x , (3.1) 



~v — 

=x 



h (Y T -Y p -C x -Cy- C c y )X p - (*;_! + k 2 )C y 



=Y 



dC e 

' ' -kiX p (El-C e x )-(k-i + h)C: 



dt 



dC 
~dt 



" -kiYp(El-C e y)-(k^ + k 2 )C e y , 



=E 2 

with initial values 



a(o) = o, <7„(o) = o, c«(o) = o, a e (o) = o. (3.2) 
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The initial values of X p and Y p are arbitrary. 



Following the approach in the previous section, we reduce Eq. (3.1) to a two dimen- 
sional system. Assuming the validity of the tQSSA, we obtain an approximating differential- 
algebraic system. Solving the algebraic equations, which are linear in the original coordi- 
nates, leads to a closed, reduced system of ODEs. We end by discussing the validity of the 
tQSSA. 



3.1. New coordinates and reduction under the tQSSA 

To extend the tQSSA we define a new set of variables by adding the concentration of 
the free state of a species to the concentrations of all intermediate complexes formed by that 
particular species as reactant [5], 



X p :— Xp + C y + C%., 
Y p := Y p + C x + Cy. 



(3.3) 



Under the tQSSA, the intermediate complexes equilibrate quickly compared to the 



variables X p and Y p . In the coordinates defined by Eq. (3.3), Eq. (3.1) takes the form 

k 2 C x — k 2 C^, 



dX n 



dt 
dYp 

dt 








k 2 Cy _ k 2 Cy, 

h(X T -X p - C X )(Y P - C x - CD - (*_! + k 2 )C x , 
h(Y T -Yp- C y )(X p -C y - C e x ) - (*_! + k 2 )C y , 
h(X p -C y - C e x )(E? - CI) - (*_! + k 2 )C e x , 
h(Y p -C x - C e y)(EZ - C e v ) - (k-i + k 2 )C e v . 



(3.4a) 

(3.4b) 

(3.4c) 
(3.4d) 
(3.4e) 
(3.4f) 



Solving the coupled system of quadratic equations (3.4c 3.4f ) in terms of X p , Y p appears to be 
possible only numerically, as it is equivalent to finding the roots of a degree 16 polynomial [3]. 
However, since we are interested in the dynamics of X p and Y p , we can proceed as in the 



Defining k„ 



(3.3 


) in ( 


3.4c 


3.4f 



Yp ~\~ k m 




Y P 







C x 




Y p (Xt — X p 


X p 


X 4- b 





Xp 




Cy 




X p (Y T — t P/ 








Xp + k m 











X V E\ 











Yp -\- k m 








Y E T 



(3.5) 



The coefficient matrix above is invertible and Eq. (3.5 ) can be solved to obtain C x , C y , C x , Cy 
as functions of X p , Y p . Denoting the resulting solutions as C X (X P , Y p ), C y (X p , Y p ), C X (X P , Y p ), 



Cy(Xp,Y p ) and using them in Eqs.(3.4a 3.4b) we obtain the closed system of equations 



d 
dt 



Xp 
Y n 



C x (X p , Yp) — C x (X p , Yp) 
Cy(X p , Y p ) — Cy(X p , Y p ) 
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Reverting to the original coordinates, X p and Y p , and using the chain rule gives 



d_ 
dt 

1 J_ 

^ dX p 
dC x , 
dX v 



X p + C y (X p , Y p ) + C*(X P , Y p ) 
Y p + C x (X p , Y p ) + Cy(X p , Y p ) 



+ 

dX v 



dC% 
dX v 



dC y 
dY v 



1 + 



+ 

dC x 
dY v 



dY v 



+ 



aci 

dY v 



d 

dt 



X p 
Y, 



C x (X p , Yp) 

Cy(X p , Yp) 

C x {X p , Yp) 

Cy(X p , Yp) 



~~ C x (X p , Y p ) 

— Cy(X p , Yp) 

~~ C x (X p , Y p ) 

- Cy(Xp,Y p ) 



. (3.6) 



The initial values of Eq. (3.6) are determined by projecting the initial values, given 



by Eq. (3.2), onto the slow manifold. Unfortunately, they can be expressed only implicitly. 



The reduction from Eq. (3.1) to Eq. (3.6) was obtained under the assumption that X v = 
X p + C y + CJ and Y p = Y p + C x + Cy are slow variables, and hence constant during the 
transient to the slow manifold. Therefore the projections of the initial conditions onto the 
slow manifold, X p (0) and Y p (0), are related to the original initial conditions as 



X p (0) + C y (X p (0), Y p (0)) + C e x (X p (0), Y p (0)) 
Y p (0) + C x (X p (0), Y p (0)) + C*{X p (0), Y P (Q)) 



X p (Q) + C y (0) + Ct{Q) 

y p (o) + c x (o) + c;(o) 



— x T , 

= Y T . 
(3-7) 

We have therefore shown that, if the tQSSA holds, and if the coefficient matrix on the 



left hand side of Eq. (3.6) is invertible, then 



d_ 
dt 



X p 
Y„ 



I _|_ ®Cy 



dX v 



+ 



dX v 



dCy 

dY v 



+ 



dC% 
dY v 



dC x 

ax v 



+ 



dX v 



1 + 



8C X 



^ dY v 



C x (X p , Yp) 

Cy(X p , Yp) 



~ C X (X P , Yp) 
— Cy(X p , Y p ) 



(3.8) 



with initial value obtained by solving Eq. (3.7), is a valid approximation of Eq. (3.1 ). Fig. 2 



shows that the solutions of the two systems are indeed close, after an initial transient. 



3.2. Validity of the tQSSA for two interacting proteins 

To reveal the asymptotic limits for which the tQSSA holds, we again rescale the original 
equations. In particular, X p and Y p are scaled by the total concentration of the respective 
proteins. To scale the intermediate complexes, each MM reaction in this network is treated 



as isolated. The scaling factors are then obtained analogously to (3 in Eq. (2.12). Let 

XjYj' 



a. 



Xx + Yx + k m 

X^pEf 
X T + Ef + k ri 



a, 



X^Y^ 
X T + Y T + k r , 

YtE t 2 
Y T + ET + K 



and 



max 



Xt Xt Yr Yr 



k 2 a x k 2 (3% k 2 a y k 2 f3y 



Therefore, T s is obtained analogously to Tx in Eq. (2.12). The reason for choosing the 



maximum will become evident shortly. The rescaled variables are now defined as 

t _ X p (t) , . Y p (t) 



:(r) :-- 



C x {t) 



a. 



Cy{T) 



X p {T) 



Cy(t) 



a, 



Xt 
cl(r) 



y P { T ) 

cm 



y t 



C £ y(t) 



(3.9) 
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Using Eq. (3.3) in the Eq. (3.1) to eliminate X P ,Y P , and then applying the rescaling 



defined by Eq. (|3.9|), to the new ODE we obtain 



ft.r 



dr 

dy P 

dr 

dc r 



k x X T Y T T s dr 



k x X T Y T T a dr 



kiX T EjT s dr 



0: 



del 



k x ElY T T s dr 



Xt 



e ^OT, 



k 2 a y T s k 2 f3 e y T t 



Y T ~° y 
[Vp- 



Y T Cyl 

hi. 



(3.10a) 
(3.10b) 



PI 



CX. X kn 



X t Yt x ^V X t Y t 



V p x T x XtV ^pyp ~ x T u xi>p "t x T vyp Yt y p i T y T S 



(3.10c) 



_ 



X T Y T V X t YtV 



Xar 



fit 

7 C~X r 



a; e 



e\2 



A.. _ 



EfXj 



C x Cy 



(3.10d) 

fi x k m 



2/p 



/3 e 



Oh 

Yi 



„e\2 



T 



EfX T 

(3.10e) 

Pyk m e 
EjYr Cy) 



where 



ko 



Y T 



h (X T + Y T + k r , 



Y2-- 



h (X T + Ef + k n 



ko 



Xn 



' y - hiYr + Xr + k^ 

* _ h El 

' y - htYT + El + kmf 



(3.10f) 



The bounds on these coefficients follow from the definition of T s . Since (1/T S ) < (A^o^/Xt), 



ft.7 



< 



k 2 at 



XjY'x 



Similarly, 



k. \ X'jYj' I $ k\ X^Y r p k\ XjYp \Xj> 4~ Yp 4~ k r 

a y < f ft* < p e j Px 

^XtYtTs-^ hXTE?T s - e *> aM k x ElY T T s 



<4 



Finally, we define 

e := m&x{e x ,e y ,e x ,e e y } . (3.11) 

The definitions of scaling factors in ( 3.9[ ) imply that all the coefficients on the right hand 
side of (3.10c-3.10f ) are Oil). Therefore, in the asymptotic limit e — > 0, Eq. (3.10) defines 
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a singularly perturbed system. Since the two equations are related by the scaling given in 
Eq. (3.9), we can conclude that in the limit e — > 0, the tQSSA is valid. If additionally the 
slow manifold is normally hyperbolic, then Eq. (3.4) is a valid reduced model of the network's 
dynamics. The normal hyperbolicity and stability of the slow manifold will be proved in a 
general setting in section [4] . 



4. The general problem 



We next describe how to obtain reduced equations describing the dynamics of a large 
class of protein interaction networks [HI EU [33l |32j El HI |25j [9] . We again assume that the 
proteins interact via MM type reactions, and that a generalization of the tQSSA holds [5]. 
We will follow the steps that lead to the reduced systems in the previous two sections: After 
describing the model and the conserved quantities, we recast the equations in terms of the 
"total" protein concentrations (cf. sections 2.1 and 3.1). Under a generalized tQSSA, these 
equations can be reduced to an algebraic-differential system. We show that the algebraic 
part of the system is linear in the original coordinates (cf. sections 2.2 and 3.1), so that 
the reduced system can be described by a differential equation with dimension equal to the 
number of interacting proteins. We next show that this reduction is justified by proving 
that the singularly perturbed system we examine satisfies the conditions of GSPT (cf. see- 
Finally, we describe the asymptotic conditions under which the system is singularly 



tion 



2.3). 



perturbed, following the arguments in sections 2.4 and 3.2 



4-1. Description of the network 

We start by defining the nodes and edges of a general protein interaction network. The 
nodes in this network represent enzymes as well as proteins, while the edges represent the 
catalytic effect one species has on another. Proteins are assumed to come in two states, 
phosphorylated and unphosphorylated. Both states are represented by a single node in this 
network. Fig. [3] and the following description make these definitions precise. 
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'E. 



(r, i) th EP-type 



edge 



r th E- typo node 



(j, i) th PP-typc 
edge 



px thPP-type ^ 

edge 





i th P- type node 




C7. 



(s,j) th EP-type 



edge 




E 



s th E- type node 



j th P- type node 



Figure 3: A simple example illustrating the terminology used in describing protein interaction networks. 
The shaded regions represent nodes and encompass either an enzyme or a single protein that is part of 
an MM type reaction. Each dotted arrow represents an edge in the network. The solid arrows represent 
transitions within the nodes, and do not define an edge in the network. 
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In a network of n interacting proteins, and n associated enzymes, we define the follow- 



ing: 



Nodes: The two types of nodes in this network represent proteins (P-type nodes) and 
enzymes (E-type nodes). Each protein can exist in either an active or inactive form. The 
inactive form of the ith protein is denoted by U, and the active form by Pj. The ith P-type 
node is formed by grouping together Ui and Pj. In addition there are n species of enzymes, 
Pj, which exist in only one state. 

Edges: All edges in the network are directed, and represent the catalytic effect of a 
species in a MM type reaction. There are two types of edges: PP-type edges connect two 
P-type nodes, while EP-type edges connect E-type nodes to P-type nodes. In particular, a 
PP-type edge from node % to node j represents the following MM type reaction in which Pj 
catalyzes the conversion of Uj to the active form Pj, 



Pi + Ui 



c 



k 2 



Pi + Pi- 



(4.1a) 



Note that autocatalysis is possible. The rate constants k\j, k^j, kfj, associated to each edge, 
can be grouped into weighted "connectivity matrices" 



K. 



IK 1 ] 



K, 



In the absence of an edge, that is, when Pj does not catalyze the phosphorylation of Uj, the 
corresponding (z,j)-th entry in K 1: K-x, and K 2 is set to zero. 

EP-type edges are similar to PP-type edges, with enzymes acting as catalysts. To 
each pair of enzyme, Pj, and protein, Pj, we associate three rate constants l]j, If j of the 
corresponding reaction in which Pj is a catalyst in the conversion of Pj into Uj, 



Ei + P 



c 



i 2 



Uj + P, 



(4.1b) 



The rate constants can again be arranged into matrices 

with zero entries again denoting the absence of interactions. 

These definitions imply that the active form of one protein always catalyzes the pro- 
duction of the active form of another protein. This assumption excludes certain interactions 
(see section [5] for an example). However, the reduction is easiest to describe under these 
assumptions, and we discuss generalizations in the Discussion. 

For notational convenience we define U = [Ux, U2, ■ ■ ■ , U n Y, P = [Pi, P 2 , 
P = [Ei,E 2 , . . . , E n Y, and arrange intermediate complexes into matrices, 



PJ*, and 



Wi 

'21 



On 



°12 

CP 



12 
•U 
22 



°ln 
U 2n 



C 



u 
111 



°n2 



C 



U 
nn 



r<E 
°n 


CE 
L/ 12 


CE 
U Xn 


°21 


r<E 

°22 . 


r<E 

°2n 


r<E 

nl 


U n2 


fiE 

nn 
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Initially all intermediate complexes are assumed to start at zero concentration. Therefore, 
any intermediate complex corresponding to a reaction that has zero rates, will remain at 
zero concentration for all time. 

For instance, in the two protein example analyzed in section [3} we have 



c 



u 



Cy 

Cy 



c F 



ct 
o ct 



u 



X 
Y 



P 



X p 

Y, 



etc. 



Assuming that the system is isolated from the environment implies that the total 
concentration of each enzyme, Ej , remains constant. Therefore, 



% 
it 



E>. + J2 C ° = E ^ «e{l,2,...,n}. 



(4.2a) 



Similarly, for each protein the total concentration, Uf , of its inactive and active form, and 
the intermediate complexes is constant, 



U i + Pi+[£c% + £c%-C%)+ y Ec* = UT, ie {1,2, ...,n} 



(4.2b) 



r=l 



r=l 



Let 



V n = [l 1 ... l] f , E T =[E T l El ... El]\ and U T = [ U( Uf ... U, 



T It 



TT TTT 



jT it 
'n J 5 



n times 



and denote the n x n identity matrix by /„. In addition, we use the Hadamard product of 
matrices, denoted by *, to simplify notation^} Constraints (4.2) can now be written concisely 
in matrix form 

E T = E + C E V n , 

Ut = U + P + Cu V n + C\jV n - (I n * Cu)V n + C%V n . 



Applying the law of mass action to the system of reactions described by (4.1a 4.1b) yields a 
(2n 2 + n) dimensional dynamical system, 



( ~~ ki S PiU s + (k is l + k1 s )C^ s j + ^ ( ferjCri — lliE r Pi + Ki^fi j , Pi 
s=l ^ ' r=l ^ ' 



dP 
dt 

dC u - 

= ^ijPiUj — {hlj + fcyjCy. 



(0) = Pi 



dC E - 

—jj- = lijEiPj - + ZyJCy, 



Cg(0) = 0. 

(4.3) 



c|(o) = 0, 



2 For instance, the Hadamard product of matrices A 

ae bf 
eg dh 



and B 



e f 
g h 



, is A * B 
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Due to the constraints (|4.2a|4.2b), Ui,Ei, are affine linear function of Pi,C^,C^ and can 







be used to close Eq. (4.3). Our aim is to reduce this 2n + n dimensional system to an n 



dimensional system involving only Pi . 

4-2. The total substrate coordinates 

In this section we generalize the change of variables to the "total" protein concentra- 



tions, introduced in Eq. (3.3). Let 



P<:=Pi + £c£ + X;C& •e{l,2,...,n} J 



(4.4) 



s=l 



r=l 



so that Eq. (4.3) takes the form 



dP\ 
dt 

dCf 3 

dt 
dCE 



'ri i 



^ j ^rfiri ^ 1 Kfir 
r=l r=l 



dt 



(4.5a) 

(4.5b) 
(4.5c) 



To close this system we use Eqs. ( 4.2a|4.2b~ ) with Eq. (4.4), to obtain 

n n 

s=l r=l 
n 

r=l 

n 

Ei = 

s=l 

n n 
Pi = Pi~ ^ ^ ~~ ^ 



(4.6) 



s=l 



r=l 



Defining P := (P 1; P 2} P n Y, Eq. (4.4) can be written in vector form as P = P + CjjV n 
C E V n , and Eqs. (4.5) and (4.6) can be written in matrix form as 



dP 

~dt 
dC v 

dt 
dC E 
~dT 



--{k 2 * CuYv n - (l 2 * c E yv n , 

-K x * {PU l ) - (K-i + K 2 ) * C v , 
* {EP l ) - + L 2 ) * C E , 



(4.7a) 
(4.7b) 
(4.7c) 
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where 



U 



E 



U T -P- C V V n - CjjVn - C E V n + (/ n * C V )V n 
■UT-P-C'uVn + iln^C^Vn, 



P — P — C\jV n — C E V n . 



(4.8a) 
(4.8b) 
(4.8c) 



4-3. The tQSSA and the resulting reduced equations 

The general form of the tQSSA states that the intermediate complexes, Cu and Ce, 



equilibrate faster than P. This assumption implies that, after a fast transient, Eq. (4.7) can 
be approximated by the differential-algebraic system 



dP 
~dt 



{k 2 * CuYv n - (l 2 * c E yv n , 

=K X * {PU l ) - + K 2 ) * Cu, 
=Li * (EP f ) - (L_i + L 2 ) * C E - 



In particular, according to GSPT (see section 2.3), if the slow manifold 

(P, Cu, Ce) 



Mo 








K x * (PU*) 
L x * (EP*) 



K 2 ) * Cu; 
L 2 ) * C E 



(4.9a) 

(4.9b) 
(4.9c) 



(4.10) 



is normally hyperbolic and stable, then the solutions of Eq. ( |4.7 ) are attracted to and shadow 
solutions onXj. 

If we consider the system ( 4.9| j,c) entry-wise then it consists of 2n 2 coupled quadratic 
equations in 2n 2 + n variables, namely the entries of P, Cu, Ce (note that U, E are functions 
of P, Cu, Ce)- As described in section 3.1, we can avoid solving coupled quadratic equations 



by seeking a solution in terms of P instead of P. Using Eq. (4.8 i,b) we eliminate E, U from 
Eqs. ( 4.9b , c) to obtain 



K x * [P (V*Cu + V l n C v - V*(I n * Cu)) + PV*C E ] + + K 2 ) * C v 

Li * (C E (VnP*)) + [L-i + L 2 ) * C E 



K x * [P (C4 - P*)] 
(4.11a) 

U * (E T P l ) . 

(4.11b) 



Although complicated, Eq. (4.11) is linear in Cu and Ce- The following Lemma, proved in 



Appendix C , shows that the equations are also solvable. 



Lemma 1. Suppose K x = [k} 3 ], K„ x = [k^], K 2 = [k%], L x = L„ x = [l^], L 2 = [Zj] 
€ ]R nxn are real matrices with non-negative entries. Furthermore, assume that for any 
pairi,j E {l,2,...,n} either kjj = fc" 1 = fcf,- = 0, or all these coefficients are positive, and 
similarly for the coefficients l}p l~ l , and tfj. IfUx, E T , P e ]R™ xl are real vectors with positive 



entries, and V n = [11 ■ ■ ■ 1]' is a vector of size n, then Eq. (4-11) has a unique solution for 
Cu,C E e R nxn m terms of P. 
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We denote the solution of Eq. (4.11) described in Lemma \U by Cu(P),Ce(P)- This 



solution can be used to close Eq. (4.9a), by using Eq. ( 4.8c[ ) to obtain 

dP dP d /- _\ d 

dt dt 



,IP + i[Cu(P)V n 



dt 

I+§p[Cu{P)V n 



+ 



+ -L(c E {P) t v n 



dP 



dP 
~dt 



(4.12) 



With Eq. (4.9a), this leads to a closed system in P, 
d 



I + Qp [Cu{P)V n 



dP 

~dt 



(K 2 * Cu{P)YV n - (L 2 * C E {P)fV n . 



(4.13) 

The initial value of Eq. (4.13), denoted by P(0), must be chosen as the projection of the 
initial value P(0) of Eq. (4.3), onto the manifold Ai . The reduction is obtained under the 
assumption that during the initial transient there has not been any significant change in 
P = P + CuV n + C E V n . Therefore the projection, P(0), of the initial conditions onto the 
slow manifold is related to the original initial conditions, U(0), -P(O), 0/(0), C E (0), by 

P(0) + Cu(P(0))V n + C|(P(0))y n = P(0) + Cv(Q)V n + C E (0)V n = P(0). 



In summary, if tQSSA is valid, then Eq. ( |4.13[ ) is a reduction of Eq. (4.3). We next 
study the stability of the slow manifold A^o defined by Eq. (4.10). This is a necessary step 



in showing that GSPT can be used to justify the validity of the reduction obtained under 
the generalized tQSSA. 



4-4- Stability of the slow manifold 

We start by introducing several definitions and some notation to simplify the compu- 



tations involved in showing that the slow manifold Aio, defined by Eq. (4.10), is normally 
hyperbolic and stable. The results also apply to the slow manifolds discussed in sections [2] 
and|3j as those are particular examples of Aio. 

Suppose that A and B are matrices of dimensions n x k and n x I, respectively. We 
denote by [A : B] the nx(k + l) matrix obtained by adjoining B to A. We use this definition 
to combine the different coefficient matrices, and let 

C:=[C V : C E ], Q 1 :=[K 1 : L[], Q 2 

We also define 



+ K 2 : L l _i + L\ 



Z :-- 



' u ' 


Z:= 


' U T -P' 


E 







jn 



In 





and 



Vo 



2n 



[1 1 ... 1]' 



2n times 



Using this notation the right hand side of Eqs. (]4.8a||4.8b|) can be written as 

u i r u T - p 

E Ej- 



z 
z 



c t v n + 

(C* 



In 





+ 



(/„ * C\j)V n 




*C t )V T 



l2n * 
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and Eq. (4.8c) can be written as P = P — CV-m- Therefore, Eqs. (4.7b 4.7c) can be merged 
to obtain 

An 

Q 1 *(PZ t )-Q 2 *C. (4.14) 



dt 



--F(C) 



The manifold A4o is defined by 

M = {C e W nx2n | Q 1 * (PZ l ) -Q 2 *C = F{C) = 0} 



To show that A^o is stable and normally hyperbolic we need to show that the Jacobian, 

dF dF 
evaluated at Ai has eigenvalues with only negative real parts. We will show that 

has eigenvalues with negative real parts everywhere, and hence at all points of A4q, a fortiori. 



vnx2n 



snx2n 



is a matrix valued function of the matrix variables C. 



The mapping F : 
dF 

Therefore —- represents differentiation with respect to a matrix. This operation is defined 
gC 

by "flattening" a rax n matrix to a mn x 1 vector and taking the gradient. More precisely, 
suppose M = [Mi : : . . . : M n ] is a m x n matrix, where M_j is the jth column of M. 
Then define 



vec (M) :-- 



M.i 
M. 2 



C 



mnx 1 



and M := diag(vec(M)) 



C 



mnxmn 



(4.15) 

Therefore, vec (M) is obtained by stacking the columns of M on top of each other, and M 
is the mn x mn diagonal matrix whose diagonal entries are given by vec (M). 

Suppose G : C pxq ->■ C mxn is a matrix valued function with X G C px<? ^ G(X) G 
C mx " Then the derivative of G with respect to X is defined as 



dG _ dvec{G) 
OX '~ d vec (X)' 



(4.16) 



where the right hand side is the Jacobian |19J . In the appendix we list some important 



properties of these operators which will be used subsequently (see Appendix B ) 



A direct application of Theorem ITT] stated in |Appendix B| yields 



OF _ dvec(F) 
llC ~ d vec (C) 



Q 



d vec (PZ*) 
dvec{C) 



Q: 



9 vec (C) 
: 9vec (C)' 



We first assume that all the entries in the connectivity matrices are positive, so that 
all entries in the matrix C are actual variables. At the end of Appendix D we show how to 
remove this assumption. 
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Replacing dvec {C)/dvec (C) with the identity matrix, J 2n 2, adding Q 2 to both side, 
using Theorems |8j 9flO TTJ and treating P and Z as independent of C we obtain 



Q 2 + 



9 vec (F) 
dvec (C) 



Qi 

Qi 
& 



oM +( / 2 ^p) W) 



a vec (C) 



, ~ t \ 9vec(CVo n ) /T „ v 



9 vec (C) 

dvec -I^C^K)* 



9 vec (C) 
dvec (y*C - V*{{I2 n f * C)) 



dvec (C) 



dvec (C) 



2n 



9 vec {C) 
] dvec (C) 



(/2 "® P) \ 9vec(C) ^e7(^ J 



Qi [- (zv 2 f „ ® i n ) - (i 2n <g) p) | (/ 2n ® v*) - (/ 2n ® v;*) (7, 

& - {ZVt ® J n ) - (I 2 n ® P^) + (/2n ® P^) (S 



2n) 



-Qi 



(zvl ® j n ) + (i 2n ® py n *) ( / 2n2 - (/£ 



Here (/ 2n )* is the matrix obtained by applying the operator defined in Eq. (4.15) to the 
transpose of P? n . 



This computation shows that the Jacobian matrix of interest has the form 
OF 



J :-- 



dC 



-Qi 



(ZV* n ® I n ) + (J 2n ® PV*) I 2n2 - {I% nJ 



(4.17) 



The following Lemma, proved in the Appendix D[ shows that this Jacobian matrix always 
has eigenvalues with negative real part. 

Lemma 2. Suppose Z E ]R^ nxl is a 2n dimensional vector with positive entries, Y e M" xl 
is an n dimensional vector with positive entries, A, V e IR 2n x2n are diagonal matrices with 
positive entries on the diagonal. Further assume that R n and R 2n are row vectors of size n 
and 2n respectively with all entries equal to 1 . Then the 2n 2 x 2n 2 matrix 



J = A 



{ZR 2n ® I n ) + (7 2n ® FP n ) J 2n2 - (J, 



+ r 



(4.18) 



has eigenvalues with strictly positive real parts. 



This Lemma applies to connectivity matrices with strictly positive entries. In | Appendix 



D.2 we show how to generalize the Lemma to the case when the connectivity matrices contain 



zero entries. In this case only the principal submatrix of the Jacobian, J, corresponding to 
the positive entries of the connectivity matrices needs to be examined. Since any principal 
submatrix of J inherits the stability properties of J, the result follows. We therefore obtain 
the following corollary. 
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Corollary 3. The manifold A4q defined in Eq. (4.10) is normally hyperbolic and stable 



4-5. Validity of tQSSA in the general setup 

We next investigate the asymptotic limits under which the tQSSA is valid in the general 
setting described at the beginning of this section. We follow the approach given in the 
previous sections to obtain a suitable rescaling of the variables. While this rescaling does 
not change the stability of the slow manifold, AAq, it allows us to more easily describe 
the asymptotic limits in which the timescales are separated, and the system is singularly 
perturbed. 



Recall that Eq. (4.7) and Eq. (4.5) are equivalent. The concise form given in Eq. (4.7) 



was useful in obtaining a reduction and checking the stability of the slow manifold. However, 



to obtain sufficient conditions for the validity of the tQSSA, we will work with Eqs. (4.5) 



and (4.6). 



Let 1% := (l^ 1 + ID/Ijj, k™ := (k^ 1 + kf^/kjj denote the MM constants. Then the 
following scaling factors are natural generalizations of those introduced in section [3j 

EjUj UfUj 
ET + Uf + lf a ' j: 1 7- I J- »,JG{l,2,...,n>. 

Note that for each pair either all of k\^ k~^ , kf^ are all zero or all nonzero. In the case 
that kjj = k'/ = jfe?. = we define k% := 0. Similarly, if Ijj = l^ 1 = 1% = then 1% := 0. Let 



l> Uj W 

T := max { max <( ) , max { ) ) = - | — , for some i ,j G {1,2, ...,n}. 



We next define the following dimensionless rescaling of the variables in Eq. (|4.5|) 

r = — , and pi{r) = — w , c£(r) = — , c?-(r) = — ^ — , !,j 6 {1, 2, n). 

J-U 1 i a ij A 



(4.19) 



After rescaling, Eqs. (4.5) take the form 



r=l \"joio Kl oJo w i "ioio^Wo^i 
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dcjj 
dr 



1-4- 



E^ e 



s=l 



E 



T is 



1 — X 



1 - 

~ 77^ E 



J r=l \ 



0!j s Cj s 



5=1 



r=l 



5=1 

1^3 



i ;/'-, + E j '/o + E + E ^ 



r=l 



5=1 



+ 



EI EJUJ 



(4.20c) 



The rescaled form of Eq. (4.5b) is similar to the rescaled form of Eq. (4.5c), and we therefore 
omit it. If we define 



€ '.i 



0? 



/ 2 



and let 



e := max < max {ej,} , max {e? } > 
^ <J hi J 



(4.21; 



then the following theorem defines the conditions under which Eq. (4.20) defines a singularly 



perturbed system and, hence, conditions under which GSPT is applicable. 

Theorem 4. If for all non-zero k\p kfj, k^ 1 and for all non zero IjjJij, Ijj and for all U?, Ej 



O 
O 

O 

O 



k 'J 

n rs 
h 2 

4.2 



o 
o 

o 

o 




o(i), 
o(i), 

0(1), 



1 < i,j,r,s < n, 



in the limit e — > ; then Eq. (4-20) is a singularly perturbed system with the structure of 
Eq. (2.8). In particular, the pi are the slow variables, and the Cij and cf- are the fast variables. 



Proof. For each % there always exist indices r, s such that fc^ ^ k 2 si . Hence, the the 



right hand side of Eq. (4.20a) is not identically zero for any i e {1,2, ...,n}. Furthermore, 



by assumption all coefficients on the right hand side of Eq. (4.20a) are 0(1) as e — > 0. This 



implies that e times the right hand side of Eq. (4.20a) is identically zero, in the limit e — > 0. 



Secondly, the definition of implies that all coefficients on the right hand side of 



Eq. (4.20c) are less than or equal to 1. Also, by definition, at least one coefficient has value 
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exactly equal to 1. Hence, the right hand side of Eq. (4.20c) is not identically zero in the 
limit e — > 0. 



The definitions of e, oiij, /3ij,Tfj imply that coefficients of in Eq. (4.20c) are less than 
or equal to e. For example 



v 



< 



l 2 B- 



e e - < e. 
y — 



Hence, in the limit e — > 0, the left hand side of Eq. (4.20c) vanishes while the right hand side 



does not. To conclude the proof we only need to show the stability of the slow manifold in 



rescaled coordinates. But we have already shown that for unsealed coordinates in section 4.4 



and a non-singular scaling of variable, as in Eq. (4.19), will not affect the eigenvalues of the 
Jacobian. □ 



Hence, under the assumptions of the above theorem, Eq. (4.20 ) has the form of Eq. 2.8a 



Hence, switching back to unsealed variables we conclude that in the limit e — > 0, tQSSA is 
valid, 



i.e. 



the reduction from Eq. (4.7) to Eq. (4.9) is valid. 



4-6. The assumption of zero initial concentrations of intermediate complexes and the choice 
of scaling 

Before concluding, we discuss the significance of zero initial concentrations of interme- 
diate complexes and the benefit of the choice of scaling we used to verify the asymptotic 
limits in which the system is singularly perturbed. Proposition [5] below proves that if the 
reaction starts with zero initial concentration of intermediate complexes then the solution 
of both Eqs. (4.7) and (4.20) are trapped in an 0(1) neighborhood of the origin. Hence, 



separation of time scale in Eq.(4.20), implied by Theorem 4 can be used to obtain the reduc- 



tion of Eq. (4.7) given by Eq. (4.9). This is important, since GSPT would not be applicable 



if the rescaling were to send 0(1) solutions of Eq. ( 4.7[ ) to solutions of Eq. (4.20) that are 
unbounded as e — )• 0. 

Proposition 5. The 2n 2 + n dimensional hypercube Q defined by 

V ■= {{Pi},{c- J },M J }\0<p l < 1, < 4 < 2, < 4 < 2, Vi, j G {l,2,..,n}} 



is invariant under the flow of Eq. (4-20) 



Proof. By the construction of the differential equations from the law of mass action, all the 
species concentration variables can take only non negative values. This together with the 
conservation constraints (4.2b) force the Pi to take values between and Uf . Therefore 
< Pi{t) < 1, Vr > 0, provided the initial conditions are chosen in Q. 



Positivity of variables also implies that c^-(r) > 0, c^(r 
Q. So we only need to show that c" (r) < 2 and 



T 



< 2. 



> if the flow starts inside 
It is sufficient to show that 
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dc u . 
«j 

dr 

But 



< 0, and 5 



< 0, or equivalently that 



dt 



< 0, and ^ 



< 0. 




=2a, ; 



n n \ 



s=l 



r=l 




r=l 



Cg=2ay 



< 0. 



- 2a y ) (l ;/ - 2a fi ) - fc#2a y ] 



Similarly we can show that Cy is decreasing when = /%. This concludes the proof. □ 

From this we conclude that the assumptions of Theorem [4] and the zero initial values 
of intermediate complexes together imply the tQSSA. 

Finally, we combine the results of section |4.4 with Theorem [4] and Proposition [5] to 
obtain the main result of this study. 

Theorem 6. If the parameters of Eq. (4-3) are such that assumptions of Theorem^ are 
satisfied and the initial values of intermediate complexes are zero, then the tQSSA holds. 
For e defined by Eq. (4-21), there exists an e such that for all < e < e , the solutions of 
Eq. (4.1) are 0(e) close to the solutions of Eq. (4-9) after an exponentially fast transient. 
Eq. ( 4-3\ ) can therefore be reduced to the n dimensional Eq. (4-13) involving only the protein 
concentrations, Pi. 



5. Discussion 



We obtained sufficient condition for the validity of tQSSA in non-isolated Michaelis- 
Menten type reactions. We therefore significantly generalized previous approaches that 
extended the MM scheme to small networks of reactions [27], and provided a theoretical 
justification of the numerical results obtained in [5]. 

We noted that the direct application of the tQSSA to equations modeling networks 
of reactions produces a reduction that contains coupled quadratic equations. However, for 
the class of networks discussed here we were able to circumvent this problem by solving 
and equivalent linear system. Moreover, we obtained a closed form equation in terms of 
protein concentrations only. A direct application of the tQSSA leads to a reduced system 
that involves the concentration of proteins and intermediate complexes. It was also shown 
that the slow manifold used in the system reduction is always attracting. 
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MM type reactions are often used in models of signaling networks. In such models it is 
frequently assumed that the reduced equation describing the dynamics of a single, isolated 
protein can be used to study interactions in networks. It has been noted that this use of 
MM differential equations is not necessarily justified [3]. The present approach provide an 
alternative approximation that was proved to be valid. 

Recently, a general reduction procedure for multiple timescale chemical reaction net- 
works has been proposed |15j . That study considered a general chemical interaction network, 
with a pre-determined set of fast and slow reactions. We deal with a more restrictive class 
of equations, which makes it unnecessary to start with a prior knowledge of fast and slow 
reactions. Moreover, we are able to show the normal hyperbolicity of the slow manifold in 
our reduction, something that was not possible in the more general setting described in [UJ . 

We end by pointing out a couple of limitations of this work. Firstly, not all enzymatic 
networks belong to the class we have considered here. For example, our full reduction scheme 
does not work for the network depicted in Fig. |4j 
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Figure 4: A hypothetical network for which the reduction described in sections |4 . 2\\4 . 3| leads to a differential- 
algebraic system of equations. The concentrations of the intermediate complexes appear in a nonlinear way 
in the resulting algebraic equations. A further reduction to a form involving only the protein concentrations 
is therefore not apparent. 
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This network is a slight modification of the network in Fig. [2^,). Although the tQSSA 
can be justified, the algebraic part of the reduced equations cannot be solved using our 
approach. These equations have the form 

= (X T — X p — C e x — C x — C y ) (Y T — Y — C y — C x — Ct) —k m C x , 
v v '> v ^ 

=x =Y P 

= X p (Y T -Y-Cy-C x -C e y )-k m Cy, 

S V ' 

=Y P 

= {El-Cl)X p -k m Cl 
= (Ej — Cy)Y — k m Cy, 

which has to be solved for C x ,C y ,C x , Cy in terms of X p , Y. Immediately we run into prob- 
lems because the first equation in the above algebraic system is quadratic in the unknown 
variables. 

We also note that no approximation theory is truly complete unless error bounds are 
investigated. Although GSPT guarantees that the derived approximations are (9(e) close to 
the true solutions, a more precise description of the error terms may be desired. 

Acknowledgements: We thank Patrick de Leenheer, Paul Smolen and Antonios 
Zagaris for helpful discussions and comments on earlier version of the manuscript. This 
work was supported by NSF Grants DMS-0604429 and DMS-0817649 and a Texas ARP /ATP 
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Appendix A. Bound on the expression for e as obtained in Eq. (2.14) 



Lemma 7. (Bound on e): If ki, &2, k-i, e, x G M+, then 

k 2 e 



e 



h ( e + x + ^±^) 2 



Proof. Since k\,k 2 , fc-i, e, x are all positive, 



< 



k\e k<i 1 
" (he + h) 2 ~ 4' 



k\e k 2 



k 1 (e + x+ k -^) 2 ~ hie + ^y (he + k 2 y 
Since for any positive number s, s + 1/s > 2, we obtain 



k\e k 2 



(he + k 2 ) 2 



< 



1 



— + 

&2 V fcie 



□ 



This bound is sharp because for k\ — 1, k 2 — 1, k-\ — > 0, e = l,x — >• we obtain 
e -> 1/4. 
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Appendix B. Differentiation with respect to a matrix 



The theory of differentiation with respect to a matrix is described in [19J. We already 
introduced the vec and hat operators and the definitions of differentiation with respect to a 



matrix variable in Eqs. (4.15) and (4.16). Below we list some important properties of these 



operators as they relate to differentiation with respect to a matrix. Proofs can be found 
in [19]. 

Theorem 8 ([281 E2])- For any three matrices A,B and C such that the matrix product 
ABC is defined, 

vec (ABC) = (C* g> A) vec (B). 

Theorem 9 ([IS]). For any two matrices A and B of equal size 

vec (A * B) = A vec (B) = Bvec (A). 

Theorem 10 (Product rule[19]). Let G : C pxq -»■ C mxr and H : C pxq -»■ C rxn be two 
differentiable function then 



d vec (GH) 



(H* 



dvec(G) 



+ (J„ ® G) 



d vec (H) 



dvec(X) v m, dvec(X) v " J dvec(X)' 

Theorem 11 (Hadamard product rule PJ). Let G : C pxq C mxn and H : C pxq 
C mxn 6e two differentiable functions then 



dvec (G*H) _ - 9vec (G) - 9vec (^) 
9vec(X) 9vec(X) + d vec (X) " 



Appendix C. Proof of Lemma [T] 



Note that the unknowns in Eq. (4.11) are matrices and the structure of the equation is 
somewhat similar to a Lyapunov equation, AX + XB = C. A standard approach to solving 
Lyapunov equations is to vectorize the matrices (see [13]), resulting in an equation of the 
type [(I m <S> A) + (£>' <8> I n )] vec (X) = vec (C). Proving solvability then essentially reduces 
to proving the non-singularity of the coefficient matrix [(I m £g> A) + (B f ® /„)]. We will use 



this approach to show the solvability of Eq. (4.11). 



In the proof of this Lemma we first assume that all possible reactions occur at nonzero 
rates so that all entries in the matrices K\, K2, if-i, Li, L2 and L_i are strictly positive. 
The result is then generalized to the case when some reaction rates are zero, so that no all 
reactions occur. 



Note that Eq. (|4.11b|) is uncoupled from Eq. (|4.11a|). Using Theorems |8| and |9| from 
section 



Appendix B[ we vectorize Eq. (4.11b) to obtain 



vec [Li * (C E (VnP*)) + (L_ 1 + L 2 )*C E ] 

= vec [L, * (C E (VnP*))] + vec [(L^ + L 2 ) * C E \ 
= U vec [C E (VnP 1 )] + + L 2 ) vec (C E ) 
= L x (PVl ® J n ) vec (C E ) + (L_! + L 2 ) vec (C E ) 
Lx (PF n * ® 4) + (£_! + L 2 )l vec (C B ) 



(c.i; 
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The following lemma shows that the matrix multiplying vec (Ce) in this equation is invert- 
ible. 

2 2 

Lemma 12. If A,B G xra are diagonal matrices with positive entries on the diagonal, 
Y G IR™ xl is a column vector with positive entries , V n — [1 1 • • • 1]* is a column vector of 
size n, and I n is the n x n identity matrix, then the n 2 x n 2 matrix 

D = A (YV* ®I n )+B 

is invertible. 



Proof. Invertibility of D is equivalent to invertibility of B 1 D. Therefore it is sufficient to 
prove the result with B = I n 2 xn 2 =: I, so that D = A (YV* <g> I n ) +1.11 A (YV* <g) I n ) does 
not have —1 as an eigenvalue, then D cannot have as an eigenvalue. Demonstrating this 
will complete the proof. Let 





A 1 




yi 


A = 


A 2 




V2 




, Y = 




A n 




Vn 



where A{ G IR" xn , i G {1, 2, ...,n} are diagonal matrices, and yi G 



yi yi yi 
V2 y-i 

Vn yn 



This implies that 



A(YV*®I n ) 

ynA n y n A 

Suppose A is an eigenvalue of A {YV^ ® I n ), and 



y x A x y x A x 
y 2 A 2 y 2 A 2 



X 



x x 
x 2 

X n 



Now 



yi l n yJn yJn 

Vn^n Vn^n 



y2A 2 
ynA n 



(C.2) 



Xi G C r , i G {1,2, ...,n} the corresponding eigenvector. Using Eq. (C.2) we have 



y x A\ y\A x y x A x 
y 2 A 2 y 2 A 2 y 2 A 2 

yn-^-n ynA n y n A n 









" X, " 




x 2 




x 2 






= A 






x n 




x n 
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This implies that for all k G {1, 2, n}, 



Ml* 


Ml* 


Mi* 




"^1* " 






M2* 




M2* 




^2 fc 


= A 










UnA^n k 












X n k 



(C.3) 



where Ai k is (k, k)-ih entry in the matrix A\, and Xj fe is the fcth entry in the vector Xj. 

Therefore, if A is an eigenvalue of A (YV£ ® J n ) then it must be an eigenvalue of one 
of its n x n principal submatrices which have the form of the coefficient matrix in Eq. (C.3) 



and whose eigenvalues we know are either zero or ^h=i Di^ik ( see reason in the footnot^J). 
Hence A can not be —1, and hence D cannot have a zero eigenvalue. □ 



This settles the problem of solvablity of Ce in Eq. (4.11b). We can use this solution to 
eliminate Ce from Eq. (4.11a). Rewriting Eq. (4.11a) with all the known terms on the right 
hand side we obtain 



K 1 *[P(V^C t u + V^C u VZ(I n *C u ))]+(K_ 1 + K 2 )*C u 

= K!*[P (C4 - P*)] - K x * [PV*C E ] . 



(C.4) 



We can write 

vec [P (V*C V + V^Cu - \%(I n * Cu))] = (I n ® P) vec + V^C V - V*(I n * C u} 

(C.5) 

Since (Cj/Vn)* is a row vector, we have vec [(C[/Ki)*] = vec (C(jV n ). Therefore, using Theo- 
rems M and [9] 

vec(\/'C^) = vec(C u V n ) = (V^I n )vec(C u ), 
veciyfcu) = (I n ®V*)vec(Cu), 
vec (V*(I n * Cu)) = (I n <g> V*) vec (I n * C v ) = (/„ <g> V# J n vec (Q,). 



Plugging these in Eq. (C.5) we get 

vec[P(K'^ + K^ - V*(In*Cv))] 

= (I n <g> P) \(V^ ® I n ) + (I n <g> ^) - (/„ <g) V;*) /J vec (C v ) 



(i n ® p)(\/* ® i n ) + (j n ® pv;*) - (In ® PK!)/n vec (C L 



We have 



y 2 A 2k y 2 A 2k 



ViAi k 
ViA 2k 



1 4 


" 1 " 




' i 




1 


n 


i 






= y iAi >> 






1 


i=i 


i 



Since the coefficient matrix in the above equation is rank one, y%Ai k is the only non-zero eigenvalue. 
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The vectorized form of the left hand side of Eq. (C.4) is 



Kt (I n <8> P){V l n ® I n ) + {I n <g> PV*) - (I n ® P^)J n + + tf_0 vec (C| 



The following Lemma shows that the matrix mutliplying vec (Cu) in this expression is in- 
vertible. 

2 2 

Lemma 13. If A,B E 1R™ xn are diagonal matrices with positive entries on the diagonal, 
Y G IR™ xl is a column vector with positive entries, V n = [1 1 • • • 1]* is a column vector of size 
n, then the n 2 x n 2 matrix 



is invertible. 



that A = I„2. Now 



(I n ®Y) (V*®I n ) 



and 



{in ® ntf) 



2/1 2/i 

2/n 2/n 



2/1 
2/n 







) + (In < 




- (/„< 


3 4j 


+ B 






lity of D is 


equivalent to invertibility of ^4. 


We can 


therefore 


assume 




2/1 








2/1 










2/1 










2/n 








2/n 










2/n 













yi 








2/1 










2/1 










2/n 








Vn 










2/n 























2/i 








2/1 










2/1 










2/n 








2/n 










Vn 



2/i 2/i 

2/n 2/n 



2/i 

2/n 



2/i 2/i 

2/n 2/n 



2/i 

2/n 
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So, 



(I n ® YV>) (I n 2 - I n ) = 



and 





(4 c 


5 y) (K* ( 


8) 4) + (/„ 




(4 ® YV$ 


4 




yi 


yi 


2/1 




2/i 










2/i 








2/n 


Vn 


2/n 




2/n 










2/n 











2/1 







2/i 


2/1 


2/i 







2/i 








2/n 







2/n 


2/n 


2/n 







2/n 



















2/i 










2/i 




2/i 


2/i 


2/i 








2/n 














2/n 


2/n 


2/n 



Clearly, its sufficient to show the invertibility of D with yi — y 2 — ... — y n — 1. We 
examine 



" 1 


1 


1 


1 










1 





' 


1 


1 


1 


1 










1 











1 





1 


1 


1 







1 








1 





1 


1 


1 







1 



















1 








1 




1 


1 


1 








1 








1 




1 


1 


1 



yi 



Vn 



yi 



yi o 



y n o 



yi 



3/1 2/1 



Vn Vn 
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Now let 

V = [ V n . . . Vi n v 21 ... v 2n •■■ v nl ... v nn ] 
be an eigenvector of D corresponding to a zero eigenvalue. We aim to show that V = 0. Let 

B = diag [6n • • • bi n &21 • • • hn ••• b nl ■■■ b nn ] . 
Then for each i,j G {1, 2, ...,n}, 



-bijVij . 



(C.6) 



r=l 



Note that the left hand side of this equation, which we denote by —A*, is independ ent o f j. 
Her 

get 



Hence, for all i,j G {1,2, ...,n} we obtain = Using this observation in Eq. (C.6) we 



n , n . 



r=l 



-Xi, Vi G {l,2,...,n}. 



This equality can be written in matrix form as 

1 + 2^s=i b ls 

1 

612 



1 

6 2 i 



1 1 + e:=i ] 



1 

bin 



l>2s 



1 

&2n 



f>nl 

1 
6n2 





" Al ' 








A n 



The coefficient matrix is diagonally dominant along the columns, and hence invertible. This 



implies that Aj = 0, and so = 0. 



□ 



Lemmas 12 and 13 together complete the proof of Lemma [T] for the case when all the 
entries in the connectivity matrices are strictly positive. This proof can be extended to 
general connectivity matrices, as stated in the Lemma [T] in the following way. 

Suppose that some of the entries in the connectivity matrix are zero. Let, 



, , , , , , , 1 1 if /A. /••,/• /••,-, arc nonzero. 

' SU<l!tlm1 '^'■■' ) - \q if jfcl = jfc-l 1 jfc?. = 



II = [h{ij)]h=i such that I L {i,j) 



1 if l}j, lij are nonzero, 



if /J. 



0. 



(C.7a) 
(C.7b) 



Hence Ik and ix are the unweighted connectivity matrices of the reaction network. The 
matrices of intermediate variables, corresponding to existing connections, now have the form 



C I U K =I K *C U C% = I L *C L . 



(C.t 
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Replacing C v with C\j K and C E with &£ in Eq. (4.11), one can easily check that the 
solution of the non-zero entries of C]f and C^f does not depend on the zero entries of Ki, 
K 2 , K_i, Li, L 2 , L_i. This observation completes the proof of Lemma [T] 



Appendix D. Stability and normal hyperbolicity of the slow manifold A^o 

Following the approach in the previous section, we first prove the result under the 
assumption that K~i,K 2 , K_i, L±, L 2 and L_i are strictly positive. At the end of this section 
we show how to generalize the proof to the case when some of the reactions do not occur. 

First we start with a preliminary lemma. 
Lemma 14. Suppose Z £ ]Ri nxl is a 2n dimensional vector with positive entries, Y £ 



inxl 



is an n dimensional vector with positive entries, and \I> = [Vfy], T = 



onx2n 



real matrices with positive entries. Let A £ C be a complex number with nonpositive real 
part. If V = [vij] £ C" x2n is a complex matrix that satisfies the following system of linear 
homogeneous equations, 

2ri 1 n 

1 \ ^ 1 \ ^ Ipij 



"3 P=1 



2n 



y. 



l ^ l ^ 

:E^ + -E 



(A - %) 



1-3 1 



1 < i < n, 
1 < 3 < n, 



(A - %) v.. 



1 < i < n, 
n + 1 < j < 2n, 



(D.la) 
(D.lb) 



s=l r=l 

then V is the zero matrix. 

Proof. Let V = [v y ] £ C nx2n satisfy Eq. ( [mj ). We will show that v {j = for all Let 

2n 



1?, 



E 



1 < i < n, 



,Ei=i%. l<j<n, 
C 3 ■ I m 



En 
i=l 



n + 1 < j < 2n. 



Then Eq. (D.l) can be written as 

1 



Vi 



Ri H Oj 



' ? (A -%)Vij, 



Setting dj 



^(A-%), we have 



1 



aijZj 



Cj , 



1 < i < n, 
1 < j < 2n, 



1 < i < n, 1 < j < 2n. 



(D.2) 



equations in the unknowns {Ri, R 2} ...,R n , Cx,C 2 , C 2n } 



By summing Eq. (|D.2|) over i and j separately we obtain the following system of linear 

(D.3a) 



2n 



2n 



n 1 1 1 

^-E — + E — 



E — 



-C~ = Rj, 1 < i < n, 



C 



i=l 



J) 



1 < j < 2n 



(D.3b) 
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Eq. (D.3) can be written in matrix form as 



Vl £-'3=1 a 1:i 






1 




1 






zian 










sr^2n 1 
£-'3 = 1 a nj 


1 




1 




ZlOr.1 






1 




1 


-1 + i V " JL 

^ zi £-ii=l an 






2/10511 


y 


.On! 






1 




1 




-1 + ^- 

Z2n 


y^n 1 




Vn 


0„2n 





"■v — 

:=A 



- 


' Ri ' 








Rn 




Ci 











(D.4) 

We next show that the the coefficient matrix, A, is invertible. This will imply that 
.Ri = Cj = 0, Vi,j. This, together with (D.2), will force to be zero and we will be done. 



To show the non-singularity of A it is sufficient to show the non singularity of the 
product of A with a non-singular diagonal matrix 

y-i 



A 



Vn 



Z2r, 



. v-^2n 1 

-vi + E i= i ^ 




1 

an 


l 




2/n + Ej=l a V 


1 


1 

a n ,2n 


i 

an 


1 

"ni 






1 

Ql,2n 


1 




i 1 
+ 2^=1 a <j2n J 



Note that X is a complex symmetric matrix (i.e. X = X 1 ). To show the non singularity of 
X, it is sufficient to show that X has no zero eigenvalue. Assume that a is an eigenvalue of 
X and u G lR 3n a corresponding eigenvector. Break X into two Hermitian matrices, 

X + X* .X -X* 



X 



+i 



2i 

:=S :=T 

where X* is the conjugate transpose of X). Then, 

a(u, u) = (Xu, u) = (Su, u) + i(Tu, u). 
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To show that a is not zero, it is sufficient to show that (Su, u) is not zero for any O^uG R 3n . 
Note that, since S, and T are Hermitian, the terms (Su,u) and (Tu,u)) are always real. 



But since X is a complex symmetric matrix, SV 



Re(Xi 



J ij 2 2 ±u ^\**-l3Ji 

where Sy, and Xij are the (2, j)-th entries of the matrices S and X respectively, and 
is the complex conjugate of the complex number X^, and Re(Xij) is the real part of Xy. 
Therefore, 



S 



-Vn + E'=i Re 



an 



a»»l 



Re 



Re-^- 

On.2* 



i?e 



O n l 



-z 1 + ELi^^ 



Recall that ay = (A — 7^). If the real part of A is nonpositive then the real parts of 
dij are negative. This implies that Re— < for all In turn, this implies that S is 
diagonally dominant, and all the eigenvalues of S are negative and real, since S is a real 
symmetric matrix. 

Therefore (Su,u) < for all u G M. 3n , and a cannot be zero. This implies that X is 
invertible, which further implies that A is invertible. So, Ri = Cj = for Eq. (D.2) 
therefore implies that V = 0. □ 



Appendix D.l. Proof of Lemma^ 

Proof. We will prove the lemma by contradiction. Let 



Then 













yi 












Z2 


, Y = 














z = 
















z 2n 




. Vn . 












Zl 




z \I n 


Zj n 


Z\l n 


ZR 2n <g> J n = 


£2 


^2 


Z2 




z 2^n 


Z2^n 


Z2^n 








®In = 










^2n 


^2n 


Z2n 




Z2nln 


Z2nln 


Z2n^n 



















2n block columns 
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'2n 



YR n — I 2n 



yi yi yi 

V2 V2 V2 
2/n l/n l/n 



YR n 



YRn 



YR ri 



-YRn 



2n block columns 



Let R ( n ] 



[1 ••• 1 1-- 

Is everywhere else. Then, 



(I 2n ® YR n ) (l 2n , - 7 2 V) 



1 ] be a row vector with a zero in the i-th. place and 



YR^ 



YR n 



YRn 



2n block columns 



Therefore, 



ZR 2n ®I n + (I 2n ® YR n ) (l 2n . - 7 2 V) 





■ Zil n 








z n I n -\- Y R n 






z n+\In 




Zn+lln + YR n . 


z n+lln 


Z>2nln 


■ Z2 n In 


Z2nln 


■ Z2 n In + YR n 



Let 



A 



A (i) 



A (2) 



A (2n) 



r(!) 



r( 2 ) 



p(2n) 



where A^, T^, fc e {1, 2, 2n} are n x n diagonal blocks of A, T respectively. Hence 

J = 





A 12 


A 21 


A22 . 
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where 



.4 



11 



+ aWybP + r« ... Zl A(D 



.4 



12 



^A (n) 

2l A« ... zxAM 



A 



21 



^n+lA (n+1) 



^n + lA(" +1 ) 



.4 



22 



z 2n A^ ... z 2n A( 2 ") 



22. 



A (2n) 



^ ^A (n+1) 

2 2 nA (2n) + A( 2 ™)ri? n + r( n+1 ) 



Let A be an eigenvalue of J, with a corresponding eigenvector 



V 



VP) 
y(2n) 



G C 2 " 2 , where fc) 



(a) n 
v i 

(fc) 
v 2 



(fc) 



G C n , G {l,2,...,2n}. 



We will show that i>j = for all I, k. By definition of eigenvalues and using the block 
structure of J we get 

^iA {1) + AWy^Vw + r«y« 

z*A( n > E?=i + A^y^Vw + r^W™) 
^ n+1 A(" +1 ) £ 2 ™i + A^+^y^y^ 1 ) + r(" +1 )v(" +1 ) 

Looking at the above equation row by row we get 

2n 

2fc A (fe) J2 y{j) + A (fc) yi?i fc) \/ (fc) + T (fc V (fe) = AV (fc) , fe G {1, 2, .., n} (D.5a) 

2n 

2 fc A (fc) ^ j) + A (fc) 7iJ n y w + T (fc) 1/ (fc) = A1/ (fc) , fe G {n + 1, 2n} (D.5b) 
i=i 



Note that Eq. (D.5) is still in matrix multiplication form. Writing it further in terms of each 



of its rows, for each k G {1,2,..., 2n} and I G {1, 2, n}, we have (For notational simplicity 
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let (A(*>) -1 ) 

2n n 

11, 1 7, 



Vi 



3=1 
2» 



h=l 



(AO 



fc e {1, ...,n}, 



(D.6a) 



(fc) 



Now, Lemma 



( fc ) \ „,( fc ) 



fc G {n + 1, 2n}. (D.6b) 



14 



applied to Eq. (D.6) immediately yields that vj k ' = for all I, fc. This 
implies that the real part of A cannot be nonpositive. This completes the proof of stability 
of J. □ 



Appendix D.2. Stability of slow manifold in the absence of some connections 



Lemma 2 show that the slow manifold defined by Eq. (4.10) is normally hyperbolic 



and stable when all entries in the connectivity matrices are positive. We next show how to 
extend the result to the case when some reactions are absent. 

Recall the definitions of the unweigh ted c onnec tivity matrices, Ik,Il, and the associ- 



ated matrices Cff and Cpf given in Eqs. (C.7) and dC.8). Let In be a n x n matrix with 



ones at the places where L\, L\, L t _ 1 are non zero and zero where L\, L\, L t _ l are zeros. Now, 
recall the definition of C in section 14.41 and define 



C =[I V I Lt ]*C 

Then, in the sense that we only need to differentiate along the coordinates corresponding to 
positive connections, one can formally write 



<9vec (C ) 
<9vec (C ) 



Ik 








(D.7) 



Replacing C with Co in the definition of F and repeating the whole process of finding 
the Jacobian of F, now with respect to Co, and using Eq. (D.7) we obtain the new Jacobian 

dvec OF (C„)) 



Jr. 



<9vec (C 



where the matrix J is the Jacobian matrix given in Eq. (4.17) 



If the connectivity matrices 
have zero entries, then l will have zero entries in the diagonal. Therefore, some eigenvalues 
of Jo will be zero. But, this does not affect the stability of slow manifold because we only 
need to look for the stability along the directions of intermediate complexes that occur in 
the reactions. That is, we only need to look at the principal submatrix of Jo corresponding 
to the positive entries in the diagonal of Jo- Let this principal submatrix be Jq. But, since 
IqJqIq = IqJIqi we see that Jq is also a principal submatrix of J. And Jq is independent of 
zero entries in the connectivity matrices. Since Lemma [2] implies that, when all the entries 
in connectivity matrices are positive, J has eigenvalues with only negative real parts, we get 
that Jq will have eigenvalues with only negative real parts. We conclude that the results 
hold even if some entries in the connectivity matrices are zero. 



Coupled enzymatic networks - tQSSA 



43 



References 

1. Bennett, M., Volfson, D., Tsimring, L., and Hasty, J. (2007) , Biophys. J. 92(10), 3501 

2. Borghans, J., de Boer, R., and Segel, L. (1996) , B. Math. Biol. 58(1), 43 

3. Briggs, G. E. and Haldane, J. B. (1925) , Biochem. J. 19(2), 338 

4. Chock, P. B. and Stadtman, E. R. (1977) , P. Natl. Acad. Sci. USA 74(7), 2766 

5. Ciliberto, A., Capuani, F., and Tyson, J. J. (2007) , PLOS Comput. Biol. 3(3) (3), e45 

6. Davidich, M. and Bornholdt, S. (2008) , J. Theor. Biol. 255(3), 269 

7. Fenichel, N. (1979) , J. Differ. Equations. 31(1), 53 

8. Frenzen, C. L. and Maini, P. K. (1988) , J. Math. Biol. 26, 689 

9. Goldbeter, A. (1991) , P. Natl. Acad. Sci. USA 88(20), 9107 

10. Goldbeter, A. and Koshland, D. E. (1981) , P. Natl. Acad. Sci. USA 78(11), 6840 

11. Hardin, H. M., Zagaris, A., Krab, K., and WesterhofT, H. V. (2009) , FEBS J. 276(19), 
5491 

12. Hek, G. (2010) , J. Math. Biol. 60(3), 347 

13. Horn, R. A. and Johnson, C. R. (1991) , Topics in matrix analysis, Chapt. 4, pp. 
268-269, Cambridge University Press 

14. Huang, C. Y. and Ferrell, J. E. (1996) , P. Natl. Acad. Sci. USA 93(19), 10078 

15. Hyeong, C. and Othmer, H. G. (2010) , J. Math. Biol. 60(3), 387 

16. Jones, C. (1995) , In Dynamical Systems, Vol. 1609 of Lecture Notes in Mathematics, 
Chapt. 2, pp. 44-118, Springer Berlin Heidelberg 

17. Kaper, T. J. (1998) , In Analyzing Multiscale Phenomena Using Singular Perturbation 
Methods: American Mathematical Society Short Course, January 5-6, 1998, Baltimore, 
Maryland (Proc. Sym. Ap.), pp. 85-132 

18. Khoo, C. F. and Hegland, M. (2008) , ANZIAM J. 50, C429 

19. Magnus, R. J. and Neudecker, H. (1985) , J. Math. Psychol. 29(4) (4), 474 

20. Michaelis, L. and Menten, M. (1913) , Biochem. Z. 

21. Murray, J. D. (2003) , Mathematical Biology II, Springer, 3rd edition 

22. Neudecker, H. (1969) , J. Am. Stat. Assoc. 64(327), 953 

23. Noethen, L. and Walcher, S. (2007) , Nonlinear Anal. -Real. 8(5), 1512 



Coupled enzymatic networks - tQSSA 44 

24. Novak, B., Pataki, Z., Ciliberto, A., and Tyson, J. J. (2001) , Chaos 11(1), 277 

25. Novak, B. and Tyson, J. J. (1993) , J. Cell. Sci. 106(4) (4), 1153 

26. Pedersen, M., Bersani, A., and Bersani, E. (2008a) , J. Math. Chem. 43(4), 1318 

27. Pedersen, M., Bersani, A., Bersani, E., and Cortese, G. (2008b) , Math. Comput. Simulat. 
79(4), 1010 

28. Roth, W. E. (1934) , Bull. Amer. Math. Soc. 40, 461 

29. Schnell, S. and Maini, P. K. (2000) , B. Math. Biol. 62(3), 483 

30. Segel, L. (1988) , B. Math. Biol. 50(6), 579 

31. Segel, L. A. and Slemrod, M. (1989) , SIAM Rev. 31(3), 446 

32. Shinar, G., Milo, R., Martinez, M. R., and Alon, U. (2007) , P. Natl. Acad. Sci. USA 
104(50), 19931 

33. Tyson, J. J., Chen, K. C, and Novak, B. (2003) , Curr. Opm. Cell. Biol. 15(2), 221 

34. Tzafriri, A. R. (2003) , B. Math. Biol. 65(6), 1111 

35. Tzafriri, A. R. and Edelman, E. R. (2004) , J. Theor. Biol. 226(3), 303 

36. Wiggins, S. (1994) , Normally hyperbolic invariant manifolds in dynamical systems. 
Springer- Verlag 

37. Zagaris, A., Kaper, H. G., and Kaper, T. J. (2004) , J. Nonlinear. Sci. 14(1), 59 



